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INTRODUCTION 

This  final  report  presents  results  from  numerical  simulations  of  two  recent 
field  experiments,  the  1987  San  Vel  Quarry  source  experiment  in  Littleton, 
Massachusetts,  and  the  1987  PACE  reflection/refraction  experiment  in 
Arizona.  The  purpose  of  this  modeling  exercise  is  to  explore  the  use  of  large- 
scale  simulation  as  an  aid  in  interpreting  such  experiments  and/or  planning 
similar  source  and  reflection/refraction  studies  in  the  future.  The  two  cases 
were  chosen  by  virtue  of  their  ongoing  interest  to  Geophysics  Laboratory  staff 
and  their  general  relevance  to  the  Test  Ban  Treaty  verification  problem. 

The  San  Vel  Quarry  simulation  is  intended  to  provide  insight  into 
conventional  quarry  blast  source  effects  in  terms  of  shot  location  and 
sequencing,  as  well  as  local  topography  and  geology.  Because  of  overall 
geometrical  and  temporal  complexity  of  the  physical  model,  conventional 
seismic  analysis  tools  are  not  practical.  For  example,  layered  half-space 
algorithms  are  inadequate  in  the  near-field  due  to  3-D  source  and  quarry 
geometries,  while  geometrical  ray  tracing  analyses  cannot  capture  the 
diffracted  and  surface  wave  mode  converted  phases.  This  leaves  discrete 
numerical  simulation  as  the  “best"  approach  for  this  study.  The  principal  issue 
addressed  is  the  practicality  of  using  large-scale  3-D  finite  element  analysis  as 
a  method  for  understanding  and  generalizing  quarry-derived  near-field  data. 

The  PACE  experiment  simulation  is  used  to  explore  the  nature  of 
scattering  inhomogeneities  in  the  lower  crust.  Data  from  extensive  world¬ 
wide  seismic  reflection  and  refraction  profiling  efforts  have  increased  our 
understanding  of  structure  and  composition  in  the  crust  and  upper  mantle. 
However,  the  structural  cause  of  ubiquitous  lower  crustal  reflections 
commonly  observed  in  extensional  regions,  e.g.,  the  Basin  and  Range  or  the 
Rhine  Graben,  remains  unexplained.  Modeling  reflection/refraction  data 
using  numerical  scattering  simulations  offers  a  practical  means  to  investigate 
this  phenomenon.  The  1987  PACE  experiment  piovidcs  an  excellent  set  of 
iclraction  and  coincident  wide-angle  reflection  data  for  this  purpose. 
Quantifying  the  scale-length  and  magnitude  of  velocity  perturbation  in  the 
lower  crust  is  important  for  establishing  petrologic  and  rheologic  constraints, 
understanding  scattering  effects  on  body  and  surface  waves,  and  separating 
intrinsic  and  scattering  attenuation. 
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CONCLUSIONS 

Regarding  the  San  Vel  Quarry  simulations,  it  is  clear  that  virtually  all  of 
the  observed  phenomenology  can  be  included,  at  least  qualitatively,  in  the 
simulation  using  a  large-scale,  explicit  finite  element  code  without  overtly 
“tuning”  the  model.  Unfortunately,  the  size  and  operational  complexity  of  th  •' 
basic  calculation  on  the  Cray-2  makes  the  requisite  suite  of  simulations 
prohibitively  expensive  in  terms  of  cpu  time  and  data  post-processing — 
certainly  well  beyond  the  resources  available  for  this  project.  The  multi¬ 
borehole  source  region  model  contributes  much  of  the  simulation’s  size  and 
complexity,  but  the  necessary  transfer  of  large  amounts  of  synthetic  data  to  a 
remote  site  and  its  processing  are  the  largest  drain  on  resources  by  a  factor  of 
five  or  more. 

The  PACE  reflection/refraction  simulations  have  clearly  demonstrated  that 
elastic  waves  are  sensitive  to  the  type  of  random  media,  when  viewed  over  a 
wide-ranges  of  incidence  angles.  The  incidence  wave  and  coda  behave 
differently  to  differences  in  the  random  media,  e.  g.,  isotropic  versus 
anisotropic  small-scale  heterogeneities.  Further  understanding  of  elastic  wave 
scattering  in  the  crust  requires  high  resolution  recordings  of  seismic  data  at 
near-vertical  to  wide-angles.  Coherency  analysis  of  both  the  direct  wave  (e.g., 
the  PmP  phase)  and  coda,  as  a  function  of  offset  and  frequency,  are  necessary 
to  constrain  possible  models  of  small-scale  velocity  heterogeneities 
Furthermore,  a  comprehensive  understand  of  crustal  scale-lengths  and 
attenuation  requires  a  better  understanding  of  the  differences  observed  m  botn 
the  back-scattered  P-  and  S-wavefields. 


THREE-DIMENSIONAL  FINITE  ELEMENT  SIMULATIONS  OF 
RIPPLE-FIRED  QUARRY  BLASTS 


Gregory  L.  Wojcik  and  David  K.  Vaughan 
Weidiinger  Associates,  Los  Altos,  California 


1.  Introduction 

This  paper  describes  some  numerical  model  simulations  of  the  San  Vel 
Quarry  experiment  in  Littleton,  Massachusetts,  performed  during  the  summer 
of  1987  by  Stump,  et  al.  (1987)  and  GL,  MIT,  and  Boston  College 
researchers.  These  simulations  are  intended  to  support  continuing  field  work 
on  the  interpretation  of  quarry  explosion  source  effects  as  they  pertain  to  Test 
Ban  Treaty  verification  issues.  The  work  has  been  performed  under  an 
ongoing  research  program  on  seismic  wave  modeling  for  experiments 
conducted  under  the  auspices  of  the  Geophysics  Laboratory,  GL  Contract 
#F  1 9628-88-C-0067. 

The  ultimate  aim  of  this  work  is  to  gain  a  better  understanding  of 
conventional  quarry  blast  source  effects  in  terms  of  shot  distribution  and 
sequencing,  as  well  as  local  topography  and  geology.  There  are  a  number  of 
approaches  available  to  us  for  understanding  this  complex  phenomenon 
including  1)  conventional  layered  half-space  interpretation,  2)  field  data 
acquisition  and  interpretation,  3)  physical  model  experiments,  and  4) 
numerical  experiments  Since  the  quarry  is  three-dimensional  with  a  very 
nonuniform  surface,,  i.e.,  a  big  hole  in  bedrock,  conventional  half-space 
methods  are  inadequate  in  the  near-field,  leaving  field  data  and  models  as  the 
best  source  of  insight. 

The  principal  issue  addressed  here  is  the  practicality  of  using  large-scale, 
3-D  finite  element  model  analysis  as  a  method  for  both  understanding  and 
gencrali/.ing  quarry-derived  near-field  data.  Although  high-resolution,  large- 
scale  2-D  models  are  readily  analyzed  on  today’s  supeicomputeis,  comparable 
3-D  models  typically  require  at  least  an  order  of  magnitude  greater  resources 
in  terms  of  memory  and  speed.  However,  the  present  quarry  problem 
requires  an  intermediate  model  size  and  level  of  complexity  that  should  be 
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practical  on  a  modem  machine  like  the  Cray-2.  On  the  other  hand,  the 
complexity  of  the  source  region,  in  terms  of  sequentially  detonated  multi¬ 
boreholes,  is  an  extremely  challenging  problem. 


2.  Experiment 

The  physical  experiment  measured  ground  acceleration  and/or  velocity  at 
ranges  of  hundreds  of  meters  to  tens  of  kilometers  from  a  series  of  quarry 
blasts  at  the  San  Vel  Quarry  in  Littleton,  Massachusetts.  The  blasts  were 
conventional  quarrying  explosions  intended  to  fracture  and  rubbleize  rock  off 
one  of  the  quarry  pit's  faces.  The  pit  itself  was  approximately  400x800x60 
feet  at  the  time  of  the  experiments.  The  typical  blast  configuration  consisted 
of  three  staggered  rows  of  explosive-filled  bore  holes  (48-72  total),  each 
approximately  58-60  feet  deep,  along  a  portion  of  the  pit's  edge.  The  holes 
were  detonated  in  succession  along  the  edge,  so-called  ripple  firing,  with  the 
detonation  sequence  (delay  time)  chosen  empirically  to  maximize  rock 
fracturing  while  minimizing  ground  motions  felt  by  nearby  homes  and 
businesses.  Although  seismic  data  were  collected  at  various  ranges,  the  set  of 
principal  interest  here  are  accelerograms  from  instruments  scattered  around 
the  pit  within  a  circle  approximately  1600  feet  in  diameter.  Locations  are 
shown  in  Fig.  1  and  an  example  of  the  shot  array  geometry  and  ideal  firing 
sequence  is  given  in  Fig.  2.  The  three-component  acceleration  records  at  the 
sites  are  given  in  Stump,  et  al.  (1987).  Three  blasts,  each  on  a  different 
section  of  the  pit  perimeter,  were  recorded  at  these  instrument  sites. 

One  purpose  of  the  experiment  was  to  examine  effects  of  ripple  firing  on 
ground  motion  spectra.  This  type  of  detonation  produces  spectral  scalloping 
that  may  be  useful  in  discriminating  small  nuclear  explosions  from  quany 
blasts.  It  was  noted,  in  both  these  as  well  as  other  experiments,  that  the  actual 
delay  times  deviated  significantly  from  those  planned — apparently  caused  by 
repeatability  problems  with  commercial  blasting  caps.  The  question  is  what 
effect  do  these  random  delay  time  errors  and  local  site  effects  have  on  ground 
motion,  and  spectral  scalloping  in  particular?  Since  this  is  difficult  to  answer 
in  the  field,  our  research  approaches  the  problem  by  means  of  numerical 
modeling.  The  objective  is  to  numerically  simulate  the  detonation  sequence  in 
a  discrete  model  of  the  quarry,  generate  a  set  of  synthetic  seismograms  at  the 
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SHOT  3 


Figure  2.  The  borehole  array  geometry  and  exact  firing  sequence  for  Shot  3. 


scaled  instrument  locations,  and  compare  synthetic  and  measured 
seismograms.  The  calculations  are  expected  to  yield  insight  into  effects  of 
shot  timing  errors  and  local  topography  and  geology  on  observed  ground 
motions. 


3.  Finite  Element  Model 

To  perform  the  simulations,  a  three-dimensional,  explicit  finite  element 
model  of  the  quarry  was  built  and  executed  on  a  Cray-2  supercomputer  at  the 
Air  Force  Weapons  Laboratory.  The  240  x  240  x  74  element  model  (=  4.6 
million  elements  including  source  region  refinement)  represented  the  quarry 
pit  and  surrounding  area  (2160  x  2160  x  656  feet)  including  all  accelerometer 
sites.  It  was  gridded  to  propagate  100  Hz  shear  (S-)waves  without  significant 
dispersion,  since  this  was  the  upper  bound  on  instrument  frequency  response. 

The  model  required  significant  new  code  development  in  order  to 
accommodate  the  large  number  of  small-diameter,  explosive  filled  bore  holes. 
Since  the  typical  hole  was  much  smaller  than  a  free-field  element,  a 
subgridding  capability  was  developed  to  grade  the  grid  size  down  from  the 
free-field  dimension,  through  an  intermediate  zone,  to  approximately  twice 
the  hole  diameter. 

Because  the  quarry  blast  is  highly  nonlinear  near  the  explosive  array — due 
to  high  pressure,  fracturing,  and  the  resulting  large  strains  and 
displacements — conventional  linear  source  modeling  techniques  could  not  be 
used.  Instead  an  “energy  pill”  source  was  implemented,  where  the  explosive 
cylinder  and  its  immediate  neighborhood  was  replaced  after  detonation  by  a 
pressurized,  outwardly  moving  region  with  one  percent  of  the  total  energy 
(1/2  kinetic  plus  1/2  potential)  as  the  explosive  products  and  included  rock. 
Nonlinear  rock  fracturing  was  mimicked  by  tension  cutoff  in  the  element's 
material  model.  A  nonlinear  cap  model  of  the  inelastic  constitutive  behavior 
of  rock  was  considered  initially  but  never  implemented  due  to  the  dominance 
of  fracturing,  rather  than  cyclic  nonlinear  processes.  Although  the  detailed 
fracture  phenomenology  is  not  simulated  by  this  model,  it  was  deemed 
adequate  for  calculating  the  seismic  pulses  radiated  by  the  explosion,  including 
shielding  by  fractured  rock  around  neighboring  boreholes.  A  two- 
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dimensional  example  of  the  model  soon  after  detonation  of  three  sheets  (i.e., 
two-dimensional  boreholes)  of  explosive  is  shown  in  Fig.  3.  This  and  other  2- 
D  models  were  used  extensively  for  development  of  the  source  and  subgrid 
coding  implemented  finally  in  3-D. 

The  principal  difficulty  encountered  with  the  model  is  late  time  numerical 
noise  contamination.  This  is  due  to  the  calculation's  relatively  long  duration  , 
about  one  second  of  simulated  time,  and  the  high  frequency  nature  of  the 
borehole  detonations  and  rock  fracturing.  Waves  reverberate  within  the 
model  many  times  over  this  period,  and  since  the  radiation  boundary  condition 
(absorbing  boundary)  on  the  outer  sides  of  the  model  is  not  perfect,  trapped 
energy  eventually  grows  to  a  significant  level.  This  problem  was  reduced  to 
the  point  that  reasonable  results  could  be  obtained  by  moving  the  bottom 
boundary  deeper  and  introducing  a  small  amount  of  viscous  damping.  The 
damping  does  not  affect  the  seismic  signals  significantly,  but  it  does  reduce  the 
ringing  substantially.  Note  that  the  damping  introduced  is  much  less  than  a 
realistic  Q  would  imply  and  is  only  introduced  for  numerical  purposes.  This 
experience  clearly  indicates  the  need  for  better  time-domain  radiation 
boundary  conditions. 


4.  Calculations 

Two  calculations  were  done  on  the  full  3-D  model.  The  first  was 
hypothetical  and  assumed  a  single  borehole  explosion  in  the  center  of  the  Shot 
3  array.  The  second  was  a  full  simulation  of  the  Shot  3  configuration  with  72 
borehole  explosions  detonated  according  to  the  original  delay  time 
specifications.  The  latter  simulation  used  about  33  hours  of  CPU  time  to 
simulate  one  second  of  model  response,  and  required  90  million  words  of 
memory.  Synthetic  vertical  velocity  seismograms  on  a  circle  surrounding  the 
Shot  3  array  and  on  lines  through  the  array,  perpendicular  and  parallel  to  the 
quarry  face  are  shown  in  Figures  4-9. 

Figures  4-6  show  the  single  shot  synthetic  seismograms.  Figure  4 
illustrates  the  influence  of  quarry  topography  on  the  azimuthal  distribution  of 
ground  motion.  Arrival  times  at  180°  to  360°  indicate  a  prominent  P-wave 
arrival  followed  by  S-waves  and  the  dominant  Rayleigh  wave.  The  P-wave  at 
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Figure  3.  A  two-dimensional  example  of  the  source  region  finite  element  model  with  grid  refinement, 
showing  a  sequence  of  snapshots  after  detonation  of  three  borehole  lines  (sheets). 
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Figure  6.  Calculated  vertical  velocity  seismograms  on  a  line 
parallel  to  the  Shot  3  quarry  face  for  a  single  shot. 
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0°  to  150°  is  a  much  weaker  arrival  due  to  shielding  by  the  quarry  faces. 
Apparently,  the  S-  and  Rayleigh  waves  are  not  shielded  except  along  azimuths 
from  130°  to  160°,  corresponding  to  the  comer  of  the  quarry.  Figures  5  and 
6  illustrate  the  decay  of  arrivals  with  distance  from  the  shot.  Figure  5  in 
particular  shows  shielding  of  the  P-wave  by  the  quarry  faces  and  a  relative 
time  delay  for  the  Rayleigh  wave  on  the  left  side  of  the  quarry.  It  also 
indicates  significantly  stronger  motion  on  the  quarry  floor  than  outside  on  the 
natural  ground  level,  probably  due  to  shielding  by  the  right  quarry  face. 

Figures  7-9  show  corresponding  synthetic  seismograms  for  the  72  ripple- 
fired  shots.  The  duration  necessary  to  capture  the  complete  sequence  (=0.8  s) 
appears  to  be  given  by  the  single  shot  duration  (=0.14  s)  plus  the  sequence 
time  (=0.7  s).  Complexity  of  these  seismograms  is  clearly  much  greater  than 
for  the  single  shot,  with  motion  dominated  by  Rayleigh  waves  from  the  long 
sequence  of  explosions.  Note  that  in  Fig.  7,  the  azimuthal  distribution  of 
amplitudes  is  similar  to  that  in  Fig.  4,  corresponding  to  the  shielding  noted 
above.  However,  the  timing  of  highest  ground  motion  at  any  azimuth  is 
directly  related  to  proximity  and  timing  of  the  shot  sequence,  e.g.,  the  strong 
arrivals  from  170°  to  200°  at  0.4  s  to  0.6  s  correspond  to  the  later  shots 
located  closest  to  these  output  points,  e.g.,  see  Fig.  2.  Similar  conclusions 
follow  from  Figures  8  and  9.  There  does  not  appear  to  be  any  directivity 
effect  due  to  shot  interference,  particularly  since  the  sequencing  would  be 
designed  to  prevent  constructive  interference. 


5.  Comparisons 

Locations  of  the  ten  accelerometer  locations  around  the  quarry  are  shown 
in  Fig.  10.  Comparisons  of  measured  and  synthetic  seismograms  at  these 
stations  are  made  in  Figures  11-15.  Figure  11  compares  arrival  times, 
durations,  and  vertical  accelerations  at  two  instrument  locations  on  the  left  and 
right  sides  of  the  quarry.  Note  that  synthetic  acceleration  was  obtained  by 
differentiating  velocity.  Significant  errors  in  both  timing,  duration,  and 
amplitudes  are  seen.  However,  the  arrival  time  comparison  is  invalid  because 
the  measured  times  are  inconsistent  with  the  instrument  distances  from  the 
shot  array.  The  amplitudes  are  in  reasonable  agreement  at  Gauge  3  but  the 
calculations  overpredict  acceleration  at  Gauge  7.  Although  truncated  by  the 
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Figure  7.  Calculated  vertical  velocity  seismograms  at  a  radius  of  280  feet  from 
the  center  of  the  Shot  3  borehole  array  for  the  72  ripple-fired  shots. 
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Figure  8.  Calculated  vertical  velocity  seismograms  on  a  line  perpen¬ 
dicular  to  the  Shot  3  quarry  face  for  the  72  ripple-fired  shots. 
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Figure  10.  Locations  of  accelerometers  around  the  quarry. 


Figure  1 1 .  Examples  of  calculated  and  measured  acceleration  histones  for  Shot  3. 


need  to  terminate  the  calculation,  the  durations  are  too  long  in  the  calculation, 
probably  due  to  the  lack  of  a  realistic  amount  of  material  damping. 

Figures  12-14  compare  synthetic  velocity  seismograms — vertical,  north, 
and  east  respectively — with  recorded  velocities  (integrated  accelerations)  at 
the  instrument  sites.  The  amplitude  ration  is  indicate  between  the  two  sets  of 
time  histories.  In  all  cases,  the  calculated  velocities  are  too  low.  This  is  not 
surprizing  considering  that  the  energy  coupling  coefficient  was  arbitrarily 
taken  as  0.01  and  could  be  significantly  higher.  Increasing  this  coefficient 
would  not  increase  all  of  the  output  stations  uniformly  due  to  the  nonlinear 
tension  cutoff  process  in  the  source  region.  The  monochromatic  appearance 
of  the  calculation  is  due  to  the  assumption  of  uniform  shot  sequencing.  A 
more  realistic  random  pattern  would  introduce  significant  interference 
between  the  arrivals. 

Figure  15  compares  synthetic  and  measured  spectra  at  the  two  sites  on 
either  side  of  the  quarry  ahown  in  Fig.  10.  It  is  clear  that  there  is  relatively 
more  high  frequency  motion  in  the  synthetics  than  in  the  data.  This  is  an 
artifact  of  the  calculation's  very  severe,  high  frequency  source  environment, 
and  inadequate  damping — both  intrinsic  and  scattering-induced — in  the  source 
region  and  on  the  travel  path. 


6.  Conclusions 

Our  conclusion  from  this  comparison  is  that  over-all  similarities  between 
calculation  and  experiment  may  exist  for  azimuthal  or  distance  variations  due 
to  gross  shielding,  but  quantitative  comparison  is  poor  in  general.  Better 
phenomenological  modeling  is  clearly  needed  in  the  source  region,  including  a 
better  estimate  of  the  source  coupling  coefficient.  This  refined  source  model 
may  require  a  subgrid  with  a  few  million  elements  alone.  In  addition, 
intrinsic  damping  should  be  included  along  with  a  weathered  layer,  i.e.,  a  high 
near-surface  velocity  gradient.  At  the  least,  a  suite  of  elastic  calculations 
should  have  preceded  those  presented  here  in  order  to  obtain  better  insight 
into  shot  sequencing  and  topographic  effects. 
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These  results  indicate  to  us  the  practicality,  albeit  limited,  of  3-D 
numerical  simulations  in  local  geology  and  topography  as  an  aid  in 
interpreting  and  perhaps  generalizing  field  data.  However,  it  is  clear  that  only 
gross  behavior  can  be  deduced  from  calculations,  i.e.,  one-to-one  comparisons 
with  field  measurements  are  virtually  impossible.  Furthermore,  the  full  suite 
of  calculations  that  would  be  required  for  an  investigation  of  shot  timing 
errors  on  spectra  and  signal  interference  are  prohibited  by  the  limited  run 
time  available  on  the  Cray-2  used  for  this  study.  Hundreds  of  hours  would  be 
necessary  for  a  more  complete  study. 

In  addition  to  studies  of  seismic  arrivals  on  the  surface,  it  would  be  of 
interest  to  examine  the  distribution  of  downgoing  energy  as  a  means  to  predict 
the  far-field  radiation  pattern.  This  can  certainly  be  done,  even  with  the 
vertically  truncated  model  used  in  this  study.  The  difficulty  is  the  limited 
capability  for  transferring  large  amounts  of  data  from  a  remote  Cray-2 
facility  to  the  analyst’s  office.  Although  possible,  in  practice,  limited 
postprocessing  budgets  currently  limit  such  studies  to  onsite  researchers  who 
have  direct  access  to  a  large  data  space. 
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Figure  12.  Comparison  of  calculated  and  measured  vertical 
seismograms  in  and  around  the  quarry  for  Shot  3. 
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Figure  1 3.  Comparison  of  calculated  and  measured  north-south 
seismograms  in  and  around  the  quarry  for  Shot  3. 
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Figure  14.  Comparison  of  calculated  and  measured  east-west 
seismograms  in  and  around  the  quarry  for  Shot  3. 
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Figure  15.  Comparison  of  calculated  and  measured 
frequency  spectra  for  Stations  3  and  7. 
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1.  Introduction 

Over  the  past  decade,  extensive  world-wide  seismic  reflection  and 
refraction  profiling  efforts  have  added  considerably  to  our  understanding  of 
structure  and  composition  in  the  crust  and  upper  mantle.  In  certain  regions 
this  work  has  increased  resolution  of  deeper  crustal  structure  to  features  on 
the  order  of  a  few  wavelengths,  although  many  questions  remain  unanswered 
concerning  the  physical  nature  of  these  structures. 

Of  particular  interest  is  the  structural  source  of  ubiquitous  lower  crustal 
reflections  commonly  observed  in  extensional  regions,  e.g.  the  Basin  and 
Range  and  Rhine  Graben.  Modeling  reflection/refraction  data  using  numerical 
scattering  simulations  offers  a  practical  means  to  investigate  this  phenomenon. 
However,  the  likelihood  of  both  lateral  and  vertical  crustal  inhomogeneity  in  a 
realistic  model  makes  such  simulation  difficult,  if  not  impossible,  using 
traditional  approaches.  This  paper  investigates  the  problem  using  stochastic 
finite  element  models  of  the  crust  and  upper  mantle. 

To  address  the  effects  of  crustal  scattering  on  refraction/wide-angle 
reflection  data,  we  model  selected  refraction  data  [ McCarthy  et  ai,  1990]  and 
coincident  reflection  data  1  Goodwin  and  McCarthy ,  1990)  recorded  acioss  the 
Buckskin-RawhiJe  Mountains  metamorphic  core  complex  in  the  southern 
Basin  and  Range.  Stochastic  models  of  lower  crustal  velocity  structure  are 
used  to  simulate  high  frequency  (>  5  Hz)  elastic  wave  propagation  in  order  to 
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determine  how  prominent  crustal  P-wave  phases  are  affected  by  small  scale 
crustal  heterogeneities.  Quantifying  the  scale-length  and  magnitude  of 
velocity  fluctuation  in  the  lower  crust  is  important  for:  1)  establishing 
petrologic  and  rheologic  constraints  based  on  models  of  admissible  velocity 
variations;  2)  understanding  the  role  of  crustal  heterogeneities  on  body  and 
surface  wave  transmission;  and  3)  determining  differences  between  intrinsic 
and  scattering  attenuation.  The  first  two  objectives  are  accomplished  here  by 
deterministic  and  stochastic  finite  element  synthetic  seismogram  modeling  of 
the  coincident  near-vertical  reflection  data  and  the  wide-angle  refraction  data. 
The  deterministic  two-dimensional  velocity  models  used  in  the  finite  element 
simulations  was  derived  from  detailed  travel  time  and  amplitude  analysis  of 
the  1987  refraction  data  [McCarthy  et  at.,  1990). 

By  way  of  background,  reflectivity  and  finite  difference  synthetic 
seismogram  techniques  have  recently  used  novel  approaches  to  simulate  lower 
crustal  reflectivity  [ Sandmeier  et  at.,  1987;  Gibson  and  Levander,  1988). 
Modeling  refraction  data  from  the  Rhine  Graben,  Sandmeier  and  Wenzel 
[1987]  simulated  lower  crustal  reflectivity  using  a  velocity  model  consisting  of 
a  series  of  random  plane  layer  thicknesses  that  alternated  between  high  and 
low  velocities.  Multiple  reverberations  within  the  finely  layered  lower  crust 
produces  the  complex  coda  seen  in  wide-angle  refraction  data. 

Gibson  and  Levander  (1988)  explained  the  same  features  by  using  a  quite 
different  velocity  model  that  contained  random  isotropic  velocity  variations  in 
the  lower  crust.  These  two  models  represent  what  can  be  considered  “end- 
member”  velocity  structures  that  may  explain  the  origins  of  lower  crustal 
reflectivity.  Importantly,  each  model,  based  on  the  dimensions  of  the  crustal 
scatterers  and  range  of  velocity  fluctuations,  lias  differing  implications  about 
the  present  state  of  the  lower  crust  and  the  processes  that  fomi  the  crust,  such 
as  the  migration  of  crustal  magma,  emplacement  of  mantle  derived  melts,  and 
lower  crustal  ductility. 

The  origins  of  body  wave  coda  observed  in  local  earthquake  recordings 
[ Frankcl  and  Clayton,  1 987 )  and  seismic  refraction  and  reflection  profiles 
[ Sandmeier  et  at.,  1987;  Gibson  and  Levander ,  1987)  have  proven  difficult  to 
explain  using  conventional  synthetic  modeling  techniques.  Crust  and  upper- 
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mantle  structure  is  quite  complex,  containing  a  large  range  of  scale  lengths 
and  velocity  fluctuations  [Wu  and  Aki,  1989]  that  are  difficult  to  accurately 
approximate.  A  variety  of  synthetic  seismogram  techniques  have  been  used  to 
model  body  wave  coda.  Ray  theoretical  methods  [ Cerveny  et  al.,  1982]  have 
been  particularly  useful  for  modeling  small-angle  scattering  in  two-  and  three- 
dimensional  structures,  where  it  is  assumed  that  the  incident  wavelength  is 
much  larger  than  the  characteristic  length  of  the  structure.  Reflectivity 
synthetic  seismograms,  limited  to  depth-dependent  variations  in  velocity,  have 
successfully  modeled  high-frequency  (>20  Hz)  elastic  wave  coda  produced  by 
a  finely  laminated  lower  crust  [Sandmeier  et  al.,  1987).  More  general  finite- 
difference  simulations  of  seismic  wave  propagation  in  two-dimensional 
random  velocity  media  [ Frankel  and  Clayton,  1984]  have  been  used  to 
constrain  the  average  scattering  properties  of  the  crust.  Finite  difference  and 
finite  element  methods,  are  being  used  more  often  in  coda  studies  because  of 
their  ability  to  model  small-scale  heterogeneities  and  wide-angle  scattering.  In 
addition,  these  techniques  are  not  restricted  to  body  wave  propagation  but  can 
also  accurately  simulate  surface  waves  in  random  media  [Hill  and  Levander, 
1984], 

2.  Basin  and  Range  Crustal  Structure 

Geophysical  surveys  and  geologic  interpretations  of  Basin  and  Range 
structure  and  tectonics  constitute  a  large  body  of  literature.  Reviews  of  the 
regional  geology  and  geophysics  of  the  Basin  and  Range  include  those  of 
Thompson  and  Burke  [1974],  Smith  [1978],  Eaton  [1979],  Speed  [1982], 
Allmendinger  et  al.  [1987],  Paktser  [1989],  Smith  et  al.  [1989],  and  Thompson 
et  al  [1989].  The  reader  is  referred  to  McCarthy  et  al.  [1990]  for  a 
description  of  the  regional  geology  and  tectonics  in  the  vicinity  of  the  1987 
PACE  experiment. 

Since  the  mid  !960’s,  extensive  refraction  profiling  of  the  Basin  and 
Range  Province  has  established  that  the  crust  is  thin  (<  32  km  thick)  with 
upper  mantle  P-wave  velocities  between  7.8  and  8.1  km  s"l.  Reinterpretation 
and  summary  of  USGS  seismic  refraction  data  [Prodchl,  1979]  suggests  a  29  to 
35  km  crust  with  a  7.9  km  s"l  upper  mantle  velocity.  Using  Nevada  Test  Site 
(NTS)  explosions  and  quarry  blasts,  Stauber  and  Boore  [1978]  found  evidence 
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for  an  anomalously  thin  crust  and  low  upper-mantle  velocity  of  7.8  km  s'l. 
Recording  independent  refraction  data,  using  NTS  explosions  and  quarry 
blasts,  Priestley  et  al.  [1982]  concluded  that  the  crust  in  northwestern  Nevada 
is  as  thin  as  20  km.  Both  Stauber  and  Boore  [1978]  and  Priestley  et  al.  [1982] 
recorded  unreversed  profiles,  hence,  their  crustal  thicknesses  and  Pn  velocities 
must  be  considered  poorly  constrained.  Refraction  results  for  northern 
Nevada  [ Benz  et  al.,  1990]  indicates  a  crustal  thickness  of  =30  km  and  upper 
mantle  velocities  of  7.9-8.0  km  s'l,-  which  argues  against  an  anomalous  crust 
in  northern  Nevada  when  viewed  on  a  Basin  wide  average.  While  details  may 
differ,  the  southern  Basin  and  Range  velocity  structure  [ McCarthy  et  al., 
1990]  is  similar  to  that  found  in  Nevada  [Benz  et  al.,  1990],  Crustal  thickness 
is  =30  km  and  upper-mantle  velocities  range  from  7.9  to  8.0  km  s' ' . 
COCORP  reflection  profiling  throughout  the  Basin  and  Range  [ Allmendinger 
et  al.,  1987;  Hague  et  al.,  1987;  Hague  et  al.,  1986]  has  also  revealed:  1)  a 
moderately  reflective  upper  crust  with  evidence  of  planar  to  listric  normal 
faults  and  asymmetric  basins;  2)  a  highly  reflective  lower  crust  marked  by 
sub-horizontal,  discontinuous  reflections;  and  3)  a  laterally  discontinuous,  but 
moderately  reflective  Moho  [Klemperer  et  al.,- 1986). 


3.  Design  of  the  1987  PACK  Seismic  Experiment 

The  1987  PACE  seismic  experiment  was  designed  to:  1)  image  the  crustal 
structure  of  the  Buckskin-Rawhide  Mountain  metamorphic  core  complex;  2) 
determine  the  crust  and  upper  mantle  structure  across  the  Colorado 
Plateau/Basin  and  Range  transition;  and  3)  record  wide-angle  refraction  and 
coincident  reflection  data  with  the  intent  of  using  the  inherent  strengths  of 
each  technique  to  improve  the  resolution  of  lower  crustal  structure. 
Achieving  these  objectives  required  deploying  a  180  km  NE-SW  trending 
refraction  profile  (Fig.  I)  crossing  the  Buckskin-Rawhide  Mountain 
metamorphic  core  complex  and  roughly  perpendicular  to  the  regional 
direction  of  extension  [Eddington  et  al.,  1986].  The  profile  consists  of  two 
140-km-long  deployments  laid  end-to-end,  each  including  120  cassette 
recorders  spaced  -750  m  apart  [Healy  et  al ,  1982],  An  in-line  400-channel, 
industry  standard  reflection  array  with  a  station  spacing  of  -40  m  was  located 
between  shotpoints  28  and  29  m  the  transition  zone  (Fig.  1).  Both  the 
refraction  and  reflection  arrays  recorded  27  shots  from  19  shotpoints. 
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Figure  1.  Location  map  of  the  1987  PACE  seismic  transect. 
Triangles  and  stars  (inset)  are  the  principal  shotpoints 
(SP)  referred  to  in  the  discussion.  The  heavy  lines  are  the 
location  of  the  coincident  reflection  arrays.  Arizona 
COCORP  profile  is  shown  as  solid  lines  in  inset. 
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Several  shotpoints  were  shot  during  both  deployments  in  order  to  increase  the 
maximum  recorded  offset  for  the  profile.  Shots  were  spaced  ~25  km  apart 
and  charge  sizes  ranged  from  400  to  2700  kg.  Details  on  the  station  locations, 
shot  size,  and  instrumentation  used  in  the  refraction  recording  are  given  in 
Larkin  et  al.  [1989],  In  the  remainder  of  this  paper,  individual  record  sections 
will  be  denoted  by  their  shotpoint  number,  e.g.,  shotpoint  21  will  be  referred 
to  as  SP21. 

Our  interpretation  of  the  1987  PACE  seismic  data  focuses  on  developing  a 
model  of  crust  and  upper-mantle  structure  based  on  finite  element  simulation 
of  selected  seismic  refraction/wide-angle  reflection  data  that  we  feel 
characterize  the  deep  crustal  structure  of  the  southern  Basin  and  Range.  Our 
discussion  complements  other  studies  of  the  1987  PACE  seismic  data  including 
the  two-dimensional  travel-time  and  ray-theoretical  synthetic-seismogram 
modeling  of  McCarthy  et  al.  [1990]  and  the  three-component  seismic 
modeling  of  Goodwin  and  McCarthy  [  1 990]. 


4.  Observed  Seismic  Refraction  and  Reflection  Data 

This  paper  utilizes  amplitude  and  travel  time  variations  observed  in  the 
refraction  and  coincident  reflection  arrays  to  infer  detailed  characteristics  of 
lower  crustal  velocity  structure.  This  study  focuses  exclusively  on  modeling 
crustal  P-waves,  and  will  primarily  address  variations  in  PntP  travel  time  and 
amplitude  in  terms  of  iower  crustal  structure  and  scattering.  Accurately 
modeling  PmP  is  particularly  important  considering  it  is  the  most  prominent 
wide-angle  phase  that  propagates  through  the  lower  crust,  and  will  be  the 
crustal  body-wave  phase  most  sensitive  to  small-scale  velocity  heterogeneities. 
Figure  2  shows  the  refraction  record  section  (SP21)  that  will  be  modeled  in 
this  study.  Each  trace  is  low-pass  filtered  to  12  Hz  and  plotted  trace 
normalized  with  a  reducing  velocity  of  6.0  km  s''.  Time-term  corrections 
were  applied  to  correct  for  systematic  travel-time  variations  due  to  velocity 
changes  near  the  recorder  \Kohler,  1988].  The  time-term  correction  reduces 
the  data  to  a  datum  coincident  with  the  elevation  of  the  source  by  applying  a 
static  travel  time  shifts  to  each  seismic  trace.  For  the  entire  profile,  time-term 
corrections  averaged  -0.08  s  and  ranged  from  -0.46  s  to  +0.25  s. 
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crust;  PiP,  the  wide-angl^,  reflection  from  the  mid-crustal  discontinuity;  PmP,  the  wide- 
angle  reflection  off  the  crust-mantle  boundary  (Moho);  and  Pn,  the  upper-mantle  head  wave. 


The  prominent  crustal  P-wave  phases  easily  recognizable  on  the  record 
section  are  the  upper  crustal  headwave  (Pg),  the  mid-crustal  reflection  (PiP), 
the  Moho  reflection  (PmP),  and  upper  mantle  headwave  (Pn).  The  Pg  phase  is 
observed  and  correlatable  as  a  first  arrivals  to  distances  greater  than  90  km. 
Large  travel  time  variations  of  the  Pg  phase  indicate  large  differences  in  near¬ 
surface  velocity  structure  along  the  profile.  On  average,  the  apparent  velocity 
of  the  Pg  phase  is  6.0  km  s'l.  The  mid-crustal  reflection  (PiP)  is  observed  as  a 
secondary  arrival  from  =30  (1.0  s  reduced  time)  to  =120  km  (0.5  s  reduced 
time).  The  PiP  closely  following  in  travel  time  the  Pg  phase  suggests  that  the 
mid-cmstal  boundary  is  relatively  shallow  (<15-18  km),  given  that  the  upper 
crustal  basement  has  a  velocity  of  6.0  km  s'*.  This  observation  is  supported 
by  the  modeling  results  of  McCarthy  el  al.  [1990]  that  place  the  mid-crustal 
boundary  at  12  km  beneath  the  center  of  the  profile.  Following  the  mid- 
crustal  reflection,  the  PmP  phase  is  the  most  prominent  secondary  arrival  and 
is  observed  from  =70  (=3.5  s  reduced  time)  to  180  km  (0.5  s  reduced  time). 
The  lack  of  pre-critical  energy  (<80  km)  may  indicate  a  transitional  crust- 
mantle  boundary  [Braile  and  Smith,  1977].  The  upper-mantle  headwave  (Pn) 
is  observed  as  a  first  arrival  starting  at  =140  km  (0.0  s  reduced  time)  and  is 
observed  to  a  distance  of  170  km  (-1 .0  s  reduced  time). 

High-quality,  single-fold  reflection  data  recorded  by  the  coincident 
reflection  array  are  shown  in  Figure  3.  The  record  section  is  a  composite 
made  from  the  recordings  of  shotpomts  28,  36,  and  29.  Amplitudes  are 
plotted  true  after  correcting  for  spherical  divergence  and  bandpass  filtering 
between  10-45  Hz.  The  figure  indicates  the  quality  of  the  data  and  some  of  the 
main  characteristics  of  the  Colorado  Plateau/Basin  and  Range  Transition. 
Based  on  the  reflection  data  and  for  purposes  of  discussion,  the  crust  is 
described  in  terms  of  four  units,  an  upper,  middle,  lower  crust,  and  upper 
mantle.  The  upper  crust  (<  4  s  two-way  travel  time  (TWTT))  displays 
prominent  reflectors  throughout,  the  most  conspicuous  being  the  “Bagdad 
Reflectors”  discussed  by  Goodwin  el  al.  [1989],  The  middle  crust,  between 
4.0-6.0  s,  is  transparent  seismically.  The  top  of  the  lower  crust  is  marked  by 
the  onset  of  reflectivity,  beginning  at  =5.8  s  and  2  km  and  increasing  to  =6.8  s 
TWTT  at  20  km.  The  entire  lower  crust  is  highly  reflectivity.  The  reflector 
“R””,  when  used  as  a  reference  mark,  indicates  roughly  3  km  of  crustal 
thickening  in  the  direction  of  the  Colorado  Plateau.  A  prominent  Moho 
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Distance  (km) 


location.  Data  corrected  for  spherical  divergence  and  bandpass  filtered 
between  10-45  Hz.  Triangles  at  bottom  of  section  separate  that  data  from 
SP28,  36,  and  29.  Horizontal  exaggeration  is  approximately  3:1  at  6.0  km  s'l 


reflection  is  missing  from  these  data,  but  some  hint  of  it  is  seen  at  =9.2  s  and 
2.5  km.  Below  =10  s,  the  upper-mantle  appears  transparent  as  reflectivity 
diminishes.  The  top  and  bottom  the  zone  of  lower  crustal  reflectivity  is  in 
agreement  with  the  results  of  McCarthy  el  al.  [1990]. 


5.  Finite  Element  Synthetic  Seismograms 

Explicit  second  or  fourth  order  finite  element/finite  difference  methods 
have  been  used  in  a  variety  of  geophysical  applications  that  range  from 
relatively  low  frequency  strong  ground  motion  modeling  [ Vidale  et  al.,  1985] 
to  high  frequency  reverse  time  migration  of  seismic  reflection  data  [Chang 
and  McMechan,  1987],  The  success  of  the  method  is  primarily  due  to  its 
ability  to  accurately  compute  full  wavefield  synthetic  seismograms  for 
arbitrarily  complex  velocity  structures  and  complex  sources.  This  ability  to 
model  complex  velocity  structure  has  found  an  application  in  the  simulation 
and  modeling  of  seismic  coda  observed  in  reflection/refraction  data  [Gibson 
and  Levander,  1988],  regional  surface  waves  [Jih  and  McLaughlin,  1988], 
earthquake  data  [Frankel  and  Clayton,  1984],  and  teleseisms  [McLaughlin, 
1986]  using  models  of  random  velocity  variation. 

In  this  study,  we  will  use  the  traditional  Galerkin  formulation  of  the  finite 
element  method  [Zienkiewicz,  1983],  a  formulation  implemented  by  Wojctk  et 
al.  [1988].  The  finite  element  synthetic  seismograms  have  not  been  used  as 
extensively  as  finite  difference  techniques,  but  has  been  particularly  successful 
in  earthquake  source  studies  [McGowan  et  al.,  1977;  Archuleta  and  Day  , 
1980;  Geller  et  al.,  1919;  Archuleta  and  Frazier  ,  1978],  Both  the  finite 
element  and  finite  difference  method  belong  to  a  general  class  of  methods 
known  as  the  method  of  weighted  residuals  [Huchner  and  Thornton,  1 982]. 
Unlike  the  finite  difference  method,  the  finite  element  method  does  not  define 
the  velocity  and  density  model  as  a  set  of  discrete  grid  points,  but  as  an 
assemblage  of  piece-wise  continuous  subdomains  or  elements,  over  which  the 
displacement  is  defined  by  an  approximating  polynomial  [ Zienkiewicz ,  1983]. 
Reduced  to  a  linear  system  of  second  order  ordinary  differential  equations,  the 
displacements  at  the  elemental  nodes  are  calculated  from  the  displacements  at 
the  previous  two  time-steps  and  the  displacements  at  the  surrounding  nodes.  A 


34 


lucid  discussion  on  the  formulation  and  error  analysis  of  the  finite  element  and 
finite  difference  techniques  can  be  found  in  Marfurt  [  1 985]. 


6.  Realization  of  a  Random  Velocity  Structure 

An  important  aspect  in  synthetic  seismogram  modeling  of  crustal 
scattering  is  the  accurate  statistical  description  of  the  random  velocity 
structure.  While  much  work  has  been  done  on  quantifying  the  size  and 
strength  of  random  scattercrs  in  the  crust  [Frankel  and  Clayton,  1984;  Gibson 
and  Levander,  1988],  only  a  limited  number  of  data  sets  have  been 
investigated.  Therefore,  fundamental  concepts  on  the  scale  length  and  strength 
of  crustal  scattering  for  different  geological  provinces  are  not  well  known. 
Thus  far,  modeling  studies  have  typically  used  either  gaussian,  exponential,  or 
self-similar  correlation  functions  [Frankel  and  Clayton,  1984;  Gibson  and 
Levander,  1988;  Jih  and  McLaughlin,  1988]  to  generate  realizations  of  the 
velocity  field. 

Figure  4a  shows  an  isotropic  random  field  that  is  typically  incorporated  as 
the  stochastic  component  of  the  velocity  model.  In  this  example,  an 
exponential  correlation  function,  with  a  correlation  length  of  200  m,  was  used 
to  generate  the  random  field.  While  isotropic  correlation  functions  have  been 
successful  in  numerical  studies  of  seismic  coda,  they  may  not  always  be 
appropriate.  This  is  particularly  true  in  the  Basin  and  Range,  where  active 
extension  since  Cenozoic  time  has  produced  a  sub-horizontal  crustal  fabric. 
Conventional  CDP  reflection  profiling  across  the  Basin  and  Range 
[Allmendinger  et  al.,-  1987;  Klemperer  et  al.,  1986;  Hague  el  al.,  1987]  and 
individual  shot  gathers  (Fig.  3)  are  in  accord  with  the  view  of  pervasive  lower 
crustal  reflectivity  that  is  dominated  by  a  sub-horizontal  reflection  pattern. 
Such  patterns  are  probably  produced  by  thermal  and  rheologic  processes  that 
have  preferentially  deformed  the  ductile  lower  crust  in  the  direction  of 
maximum  extension.  While  conventional  CMP  processing  of  deep  crustal 
reflections  may  produce  the  appearance  of  sub-horizontal  reflectivity  [Gibson 
and  Levander,  1988],  it  is  difficult  to  argue  that  sub-horizontal  lower  crust 
fabric  do  not  significantly  contribute,  given  the  recent  tectonic  evolution  of 
the  Basin  and  Range. 
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EXPONENTIAL  RANDOM  MEDIA 

DISTANCE  (KM) 


ISOTROPIC  200  m 


ANISOTROPIC  1  600x200m 


Figure  4.  Two-dimensional  realizations  ul  (a)  an  exponential  isotropic 
random  media  with  a  horizontal  and  vertical  correlation  length 
of  200  m  and  (b)  an  exponential  anisotropic  random  media  with 
a  horizontal  and  vertical  correlation  length  of  1600  m  and  200 
m,  respectively.  The  light  and  dark  regions  represent  high  and 
low  velocities  about  a  mean  background  velocity. 
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An  alternative  model  to  the  isotropic  random  media  is  a  realization  using 
an  anisotropic  correlation  function,  where  the  spatial  lag  differs  between  the 
horizontal  and  vertical  directions.  Shown  in  Figure  4b  is  an  anisotropic 
random  media  where  the  horizontal  and  vertical  correlation  lengths  are  1600 
m  and  200  m,  respectively.  We  feel  this  model  represents  a  compromise 
between  the  velocity  structure  used  by  Gibson  and  Levander  [1989]  to  model 
the  Rhine  graben  refraction  data  and  the  one-dimensional  plane  layer  model  of 
Sandmeier  et  al.  [1987].  It  also  approximates  a  sub-horizontal  fabric  that 
might  be  representative,  in  a  generic  way,  of  the  lower  crust  in  the  Basin  and 
Range. 


7.  Finite  Element  Simulations  of  the  1987  PACE  Seismic  Data 
Wide-angle  Refraction  Data  from  SP21 

To  investigate  the  effects  that  different  random  media  have  on  wave 
propagation,  we  will  compared  observed  seismic  refraction  data  (SP21)  with 
two-dimensional  synthetic  seismogram  simulations  that  incorporate  random 
velocity  variation  in  the  lower  crust.  The  deterministic  two-dimensional 
velocity  structure  used  as  a  basis  was  derived  from  iterative  travel  time 
modeling  [ McCarthy  et  al.,  1990].  The  three  synthetic  seismogram  record 
sections  presented  were  calculated  from:  1)  a  finite  element  grid  derived  from 
the  velocity  model  of  McCarthy  et  al.  [1990];  2)  a  lower  crust  velocity  model 
with  isotropic  small-scale  heterogeneities,  similar  to  that  in  Figure  4a;  and  3)  a 
lower  crustal  velocity  model  with  anisotropic  small-scale  heterogeneities, 
similar  to  that  in  Figure  4b, 

Figure  5  shows  the  two-dimensional  velocity  structure  [McCarthy  et  al., 
1990]  used  in  the  three  finite  element  simulations.  The  two-dimensional  finite 
element  velocity  model  was  constructed  by  digitizing  the  raytrace  model  (Fig. 
5)  every  43  m.  To  avoid  prohibitively  long  run  times  due  to  low  sediment 
velocities,  P-wave  velocities  less  than  6.0  km  s'l  were  reset  to  a  velocity  of 
6.0  km  s'l.  This  enabled  the  calculation  of  high  frequency  synthetic 
seismograms,  given  a  43  m  nodal  spacing,  but  eliminated  the  effects  of  near¬ 
surface  sedimentary  structure  in  the  calculations.  S-wave  velocities  were 
scaled  relative  to  the  P-wave  velocities  assuming  a  Poisson's  ratio  of  0.25  and 


37 


38 


densities  were  calculated  using  the  empirical  relationship  p  =  0.2!  2+0.3977a 
[Nafe  and  Drake,  1957],  where  p  is  the  density  and  a  is  the  P-wave  velocity. 

The  left  side  of  the  finite  element  grid  is  located  at  SP21  and  is  modeled  as 
a  plane  of  symmetry.  The  bottom  and  right  edges  of  the  grid  are  A2 
absorbing  boundaries  [Clayton  and  Engquist ,  1977]  and  the  top  of  the  grid  is 
a  stress-free  boundary.  The  finite  element  model  was  180  km  in  length  and  37 
km  in  depth.  The  source  is  an  isotropic  line  source  and  the  source-time 
function  is  a  Ricker  wavelet  with  a  half-power  frequency  of  8.0  Hz,  resulting 
in  approximately  17  nodes  per  wavelength  for  the  slowest  P-wave  and  10 
nodes  per  wavelength  for  the  slowest  S-wave.  The  finite  element  formulation 
assumes  a  purely  elastic  model,  hence,  intrinsic  attenuation  is  not  accounted 
for  in  the  calculations. 

Figure  6  shows  a  comparison  of  the  observed  refraction  section  and  the 
finite  element  synthetic  seismograms  for  a  deterministic  medium.  The 
comparison  shows  that  the  relative  travel  time  and  amplitude  of  the  major 
crustal  P-wave  phases  are  well  modeled  by  the  finite  element  synthetic 
seismograms.  The  existence  of  significant  pre-critical  PmP  energy  in  the 
synthetic  seismograms  suggests  that  the  model's  crust-mantle  transition  is  too 
sharp  and  that  a  transitional  Moho  is  more  appropriate  [Braile  and  Smith, 
1977].  For  offsets  greater  than  60  km,  Pg  is  over-predicted  relative  to  the 
observed  data.  In  general,  the  Basin  and  Range  exhibits  low  upper  crustal  Q 
(<200-300;  Braile,  1977;  Benz  el  ai,  1990),  which  causes  noticeable  amplitude 
and  frequency  loss  for  offsets  greater  than  60  km.  Since  the  finite  element 
simulations  assume  purely  elastic  media,  intrinsic  attenuation  cannot  be 
accounted  for  and  the  resulting  synthetics  over-predict  the  Pg  and  PiP 
amplitudes  at  larger  offsets. 

Isotropic  Random  Medium — In  this  simulation,  an  isotropic,  exponential 
random  velocity  structure  was  incorporated  into  the  lower  crust.  The  random 
medium  is  similar  to  that  seen  in  Figure  4a.  The  random  medium  was 
generated  with  a  correlation  length  of  200  m  and  a  5%  standard  deviation 
(std)  relative  to  an  average  lower  crustal  background  velocity  of  6.5  km  s‘‘. 
Shown  in  Figure  7  is  a  comparison  of  the  observed  and  theoretical 
seismograms.  It  can  be  seen  that  the  scattering  effects  from  this  model  are 
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Figure  6.  Comparison  of  trace  normalized  refraction  data  (SP21)  and  finite  element  synthetic  seismograms 
computed  for  the  two-dimensional  velocity  model  shown  in  Figure  4  but  without  scatterers. 


a)  anisotropic  random  velocities  in  the  lower  crust,  b)  isotropic  random  velocities  in  the  lower 
crust.  The  random  medium  has  a  5%  standard  deviation  (std)  about  a  mean  lower  crustal  velocity 
of  6.5  km  s*l.  For  the  anisotropic  case  the  horizontal  and  vertical  correlation  lengths  are  1600 
and  200  m,  respectively,  and  200  m  for  the  isotropic  case. 


noticeable  as  coda  following  the  PiP  and  PmP  phases  The  largest  difference 
between  this  simulation  and  the  previous  one  is  the  loss  of  correlatable  pre- 
critical  PmP  energy  for  distances  less  than  60  km.  This  simulation  suggests 
that  a  transitional  Moho  is  not  required  to  explain  the  lack  of  pre-critical  PmP 
energy  in  the  observed  data.  Results  show  that  lower  crustal  velocity 
fluctuations  distort  the  amplitude  and  phase  of  the  relatively  weak,  precritical 
PmP,  giving  it  the  appearance  of  a  transitional  Moho.  Beyond  the  critical 
point,  the  amplitude  of  the  PmP  increases  such  that  coherency  is  maintained 
with  offset  but  noticeable  amplitude  decay  occurs.  This  is  due  to  significant  P- 
to-P  and  P-to-S  conversions  that  effectively  partition  energy  away  from  the 
direct  wave,  in  this  case  the  PmP  phase.  Possibly,  the  standard  deviation  of 
the  random  velocity  structure  is  too  large  and  the  consequence  is  significant 
amplitude  loss  of  PmP. 

Anisotropic  Random  Media — In  this  simulation,  an  anisotropic, 
exponential,  random  velocity  structure,  similar  to  that  shown  in  Figure  4b,  is 
incorporated  into  the  lower  crust.  The  horizontal  and  vertical  correlation 
lengths  are  1600  m  and  200  m,  respectively.  Like  the  previous  simulation,  the 
random  fluctuations  in  velocity  have  a  5%  std  relative  to  an  average  lower 
crustal  background  velocity  of  6.5  km  s'1.  The  observed  and  theoretical 
seismograms  for  this  simulation  are  shown  in  Figure  7.  When  compared  to 
the  previous  simulation,  the  synthetic  seismograms  in  Figure  7  show  quite  a 
difference  in  both  the  character  of  the  coda  and  in  the  amplitude  variations  of 
the  PiP  and  PmP  phase.  The  coda  following  PiP  and  PmP  display  significant 
coherency  that  appears  as  prominent  phases  that  can  be  correlated  up  to  20 
km.  The  best  example  of  this  is  the  phase  that  follows  closely  in  time  PiP 
between  60  and  90  km.  These  coherent  and  relative  large  coda  waves  have  the 
effect  of  decreasing  the  coherency  of  the  principal  body  waves,  both  PiP  and 
PmP.  Similar  to  the  previous  simulation,  the  precritical  PmP  reflection  is 
weak  and  incoherent. 

When  compared  to  the  observed  data,  the  synthetic  seismograms  replicate 
some  of  the  complexity  of  the  data.  The  majority  of  coda  is  seen  as  wide- 
angle  P-to-P  conversions,  which  produced  the  relatively  large  coherent  coda. 
Correlatable  coda  energy  is  observed  as  phases  that  cannot  by  described  and 
modeled  in  conventional  terms,  like  Pg,  PiP,  and  PmP.  For  example,  between 
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80-100  km,  weak  coherent  energy  precedes  PmP  by  -0.5-1. 0  s  reduced  time. 
In  addition,  an  enechelon  pattern  of  energy  is  observed  to  follow  the  PmP  for 
=2.0  s,  which  is  similar  to  that  found  in  the  synthetic  seismograms.  Unlike  the 
observed  data,  the  synthetic  seismograms  predict  a  PmP  phase  that  is  not  as 
coherent  between  80-180  km  as  that  observed.  This  observation  demonstrates 
the  sensitivity  of  the  incident  wave  to  the  correlation  length  of  the  scattering 
media  and  possibly  the  size  of  the  velocity  fluctuations.  Perhaps,  had  the 
aspect  ratio  or  standard  deviation  been  smaller,  the  simulations  may  have 
predicted  a  higher  degree  of  coherency  for  the  PmP  phase. 

These  simulations  demonstrate  that  certain  aspects  of  each  model  fit  the 
observed  data.  The  isotropic  random  media  predicts  more  coherency  of  PmP, 
but  less  coherency  of  the  coda.  By  way  of  contrast,  the  anisotropic  random 
media  predict  coherency  in  the  coda  and  less  coherency  of  PmP.  Both  these 
results  can  be  explained  in  terms  of  the  types  of  body  waves  scattered  from  the 
incident  wavefield.  In  an  isotropic  random  medium,  for  all  angles  of 
incidence,  the  incident  wavefield  will  produce  scattered  P  and  S  energy. 
Because  the  dimensions  of  the  scattering  body  are  close  to  that  of  the  incident 
wavefield,  scattering  occurs  over  a  wide  range  of  angles  and,  therefore,  no 
coherent  wide-angle  coda  is  produced.  For  the  anisotropic  random  media,  the 
strongest  scattering  will  occur  over  a  restricted  range  of  incident  angles, 
which  correspond  to  wide-angle  P-to-P  conversion.  This  process  produces  the 
coherent  wide-angle  P-wave  coda  that  has  an  apparent  velocity  similar  to 
PmP. 

These  results  represent  only  a  limited  number  of  simulations,  so  only 
general  comments  can  be  made  about  lower  crustal  structure.  This  synthetic 
seismograms  for  the  isotropic  case  suggest  that  5%  std  is  too  large,  since  PmP 
cannot  propagate  to  large  distances  due  to  effective  P-to-P  and  P-to-S 
scattering,  which  rapidly  partitions  energy  away  from  PmP.  The  anisotropic 
case  also  indicates  that  5%  std  is  possibly  too  large,  because  the  amplitudes  of 
the  wide-angle  P-wave  coda  are  over  predicted,  compared  to  PmP.  It  is  also 
possible  that  the  aspect  ratio  for  the  anisotropic  case  is  too  large,  considering 
that  the  synthetic  seismograms  predict  coherent  P-wave  coda  over  distances 
ranges  of  30  km,  which  is  larger  than  that  seen  in  the  data.  At  this  point,  not 
enough  simulations  have  been  done  for  a  range  of  velocity  models  to 
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determine  what  is  the  range  of  permissible  structures  that  predict  the  observed 
data. 


8.  Near-Vertical  Reflection  Data  from  SP28 

Near-vertical  synthetic  seismograms  have  been  calculated  to  understand 
the  near-vertical  response  of  the  two  random  media  investigated  in  the 
previous  section.  The  simulations  are  computed  using  a  velocity  structure 
constructed  by  digitized  the  two-dimensional  velocity  (Fig.  5)  every  5  m 
starting  near  SP28,  The  finite  element  model  has  a  length  of  20  km  and  depth 
of  40  km.  Peak  frequency  of  the  Ricker  source  wavelet  is  20  Hz.  Based  on 
the  calculations  in  the  previous  section,  it  is  assumed  that  a  5%  std  is  too  large, 
and  the  following  simulations  assume  a  3%  std  about  a  mean  lower  crustal 
velocity  of  6.5  km  s'l.  Figure  8  shows  examples  of  the  velocity-depth 
function  selected  from  one  location  (SP28).  Figure  8a  is  the  velocity-depth 
function  assuming  no  random  distribution  of  velocity  in  the  lower  crust,  while 
Figure  8b  is  the  velocity-depth  function  assuming  random  variations  of 
velocity  in  the  lower  crust.  In  addition,  synthetic  seismograms  are  also 
calculated  for  a  model  with  a  transitional  crust-mantle  boundary  (Fig.  8c).  A 
transitional  crust-mantle  boundary  is  investigated  since  the  lower-frequency 
wide-angle  synthetic  seismograms  (Fig.  7  )  argued  against  a  transitional  Moho, 
based  on  the  lack  of  coherent  pre-critical  PmP  energy.  These  same  wide- 
angle  synthetic  seismograms  contradicted  the  results  from  the  model  with 
smoothly-varying  lower  crustal  velocities  (Fig.  6)  and,  therefore,  proved 
inconclusive  about  the  transitional  nature  of  the  Moho. 

The  near-vertical  synthetic  seismograms  for  the  different  cases  is  shown  in 
Figure  9.  The  single-fold  shot  gathers  are  corrected  for  spherical  divergence 
and  plotted  relative  true  amplitude.  The  synthetic  seismograms  for  the 
isotropic  and  anisotropic  random  media  and  1st  order  Moho  show  that  the 
most  prominent  phase  is  the  PmP  reflection  (=9.5  s  TWTT).  As  expected,  the 
anisotropic  case  predicts  a  greater  degree  of  sub-horizontal  reflectivity.  For 
the  isotropic  case,  the  coda  is  longer  in  duration  and  is  seen  as  the  reflectivity 
arriving  after  the  PmP  reflection  (>9.5  s  TWTT).  Over  a  wide-range  of 
incidence  angles,  isotropic  random  media  will  produce  both  P-to-P  and  P-to-S 
scattering,  relative  to  an  incident  P-wave.  Since  the  travel  time  of  back- 
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a)  the  velocity-depth  structure  of  McCarthy  et  al.  [1990], 

b)  a  random  velocity  structure  in  the  lower  crust  and  1st 
order  Moho,  and  c)  a  random  velocity  structure  and 
transitional  Moho.  The  velocity  perturbation  is  3%  about  a 
mean  background  velocity  of  6.5  km  s*l. 


ISOTROPIC 


ANISOTROPIC 


DISTANCE  (km) 


DISTANCE  (km) 


Figure  9.  Finite  element  synthetic  seismograms  for  SP28.  The  two- 
dimensional  velocity  model  included  isotropic  and  anisotropic 
random  media  in  the  lower  crust.  The  random  media  is  3%  std 
about  a  mean  lower  crustal  velocity  of  6.5  km  s-*  Source-time 
function  is  a  Ricker  wavelet  with  a  peak  frequency  of  12  Hz. 
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scattered  S-waves  are  longer,  the  corresponding  coda  will  be  longer.  This 
contrasts  with  the  anisotropic  case,  which  is  dominated  by  P-to-P  scattering, 
due  to  the  steep  incidence  angles  and  sub-horizontal  fabric.  The  coda  is 
mostly  back-scattered  P-wave  energy  with  travel  times  closer  to  the  incident 
P-wave  and,  subsequently,  shorter  coda  durations.  Synthetic  seismograms  for 
the  transitional  Moho  models  are  similar  to  the  previous  simulations  except  of 
the  lack  of  a  prominent  PmP  reflection,  which  is  in  agreement  with  the 
observed  data  (Fig  3).  Based  on  a  qualitative  comparisons,  the  synthetic 
seismograms  for  the  anisotropic  random  media  and  transitional  Moho  appear 
to  best  fit  the  observed  near-vertical  reflection  data  (Fig.  3). 


9.  Discussion  and  Conclusions 

Composition  From  Seismic  Velocity  Models — The  P-wave  crustal 
structure  derived  here  provides  constraints  on  the  composition  and  structure 
of  the  upper  lithosphere.  To  place  constraints  on  crustal  composition  in  the 
southern  Basin  and  Range,  we  show  plausible  rock  types,  taken  from  the 
laboratory  measurements  of  seismic  velocities  by  Christensen  [1979]  and 
plotted  relative  to  the  generalized  velocity-depth  function  (Fig.  10A).  The 
laboratory  determined  velocities  are  corrected  for  temperature  assuming  a 
regional  surficial  heat  flow  of  90  mW  m'2  and  pressure  effects  assuming  a 
mean  crustal  density  of  2.8  g  crn‘3,  The  crustal  geotherm  used  to  correct  the 
laboratory  velocities  is  shown  in  Figure  1  IB. 

The  important  observation  (Fig.  10A)  is  that  no  one  rock  type  follows  the 
interpreted  velocity-depth  model  for  more  than  a  few  kilometers.  We  have 
chosen  only  a  small  number  of  rock  compositions  and  metamorphic  grades 
from  Christensen  [1979],  but  we  feel  these  reflect  a  representative  sampling 
of  plausible  rock  types.  We  assume  that  granulite-grade  metamorphism  is 
achieved  in  the  lower  crust  of  the  Basin  and  Range  due  to  the  high  pressures 
(900  GPa)  and  temperatures  (>900°C).  Based  on  these  observations,  the  upper 
crust  appears  to  best  fit  a  general  suite  of  granitic-dioritic  to  quartz-rich 
granulitic  rocks  (5.5  to  6.2  km  s-'),  while  amphibolite  to  mafic  granulites  fit 
the  lower  crust,  6.4  to  6.8  km  s'l.  The  correlations  imply  that  the  bulk 
crustal  composition  varies  with  depth,  but  the  correlation  of  velocity  with 
composition  is  non-unique.  The  primary  difficulty  arises  from  the  fact  that 
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grades  are  superimposed  on  the  generalized  velocity-depth  function.  The  compositions, 
metamorphic  grades,  temperature,  and  pressure  coefficients  were  taken  from  Christensen  [1979]. 
Temperature  and  pressure  corrections  with  depth  were  made  assuming  a  90  mW  m-2  surface  heat 
flow  and  a  mean  crustal  density  of  2.67  g  cm'3. 


the  lower  crust  has  a  seismic  velocity  gradient  that  is  difficult  to  constrain,  but 
likely  ranges  from  0.02  to  0.04  s'*,  whereas  most  rock  samples  exhibit  weak 
pressure  derivatives  (0.002  s-*  or  less;  Christensen,  1979)  at  pressures 
appropriate  for  the  lower  crust.  Furthermore,  the  high  lower-crustal 
temperatures  in  this  region  of  high  heat  flow  should  further  decrease  the 
velocity  gradient  at  depth  [Christensen,  1979).  A  3%  std  in  velocity  for  the 
lower  crust  tend  to  bracket  the  range  of  crustal  compositions,  implying  that 
the  lower  crust  may  consist  of  a  large  compositional  suite.  This  observation  is 
in  agreement  with  laboratory  studies  of  exposed  portions  of  the  lower  crust 
that  show  large  variations  in  velocity  and  composition  [Fountain,  1976). 

Wide-angle  Elastic  Wave  Propagation  in  Random  Media — These 
simulations  have  clearly  demonstrated  that  elastic  waves  are  sensitive  to  the 
type  of  random  media,  when  viewed  over  a  wide-ranges  of  incidence  angles. 
The  incidence  wave  and  coda  behave  differently  to  differences  in  the  random 
media,  e.  g.  isotropic  versus  anisotropic  small-scale  heterogeneities.  Further 
understanding  of  elastic  wave  scattering  in  the  crust  requires  high  resolution 
recordings  of  seismic  data  at  near-vertical  to  wide-angles.  Coherency  analysis 
of  both  the  direct  wave  (e.g.  the  PmP  phase)  and  coda,  as  a  function  of  offset 
and  frequency,  are  necessary  to  constrain  possible  models  of  small-scale 
velocity  heterogeneities.  Furthermore,  a  comprehensive  understand  of  crustal 
scale-lengths  and  attenuation  requires  a  better  understanding  of  the  differences 
observed  in  both  the  back-scattered  P-  and  S-wavefields.  Investigations  of 
crustal  scattering  have  important  implications  over  a  broad  range  of  seismic 
applications,  these  include:  1)  understanding  the  generation  and  propagation  of 
regional  phases  used  in  the  location,  discrimination,  and  yield  estimation  of 
nuclear  explosions;  2)  determining  the  bulk  properties  of  the  crust,  based  on 
composition  and  velocity  relationships;  3)  understanding  and  differentiating 
the  role  of  scattering  and  intrinsic  attenuation  in  the  crust,  which  has 
implications  concerning  the  crust's  physical  state. 
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